{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "81c89236",
   "metadata": {},
   "outputs": [],
   "source": [
    "#############################################################################\n",
    "# 正射校正测试\n",
    "# Copyright (C) 2024  Xu Ruijun\n",
    "#\n",
    "# This program is free software: you can redistribute it and/or modify\n",
    "# it under the terms of the GNU General Public License as published by\n",
    "# the Free Software Foundation, either version 3 of the License, or\n",
    "# (at your option) any later version.\n",
    "#\n",
    "# This program is distributed in the hope that it will be useful,\n",
    "# but WITHOUT ANY WARRANTY; without even the implied warranty of\n",
    "# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the\n",
    "# GNU General Public License for more details.\n",
    "#\n",
    "# You should have received a copy of the GNU General Public License\n",
    "# along with this program.  If not, see <https://www.gnu.org/licenses/>.\n",
    "#############################################################################"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "5d373774",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from osgeo import gdal, osr\n",
    "from myraster import MyRaster\n",
    "from sklearn import linear_model\n",
    "from sklearn.preprocessing import PolynomialFeatures\n",
    "from scipy.optimize import curve_fit\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "ad35c266",
   "metadata": {},
   "outputs": [],
   "source": [
    "def read_points(path):\n",
    "    with open(path, 'r') as f:\n",
    "        res = []\n",
    "        crs = f.readline()\n",
    "        assert crs[:6] == '#CRS: '\n",
    "        osra = osr.SpatialReference(crs[6:])\n",
    "        assert f.readline() == 'mapX,mapY,sourceX,sourceY,enable,dX,dY,residual\\n'\n",
    "        while True:\n",
    "            line = f.readline()\n",
    "            if len(line) == 0:\n",
    "                break\n",
    "            res.append(tuple(map(float, line.split(','))))\n",
    "    return osra, res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "5060b867",
   "metadata": {},
   "outputs": [],
   "source": [
    "# DEM是从USGS EarthExplorer下载的多张30m SRTM, 并用gdal合并的虚拟图层\n",
    "dem=MyRaster(gdal.Open('/mnt/2thdd-btrfs/rs/DEM/SRTM1_v3/USGS_EE/virtal.vrt'))\n",
    "# 地面控制点文件是用QGIS的Georeferencer工具制作\n",
    "# USGS解密影像与现代资料（OpenStreetMap，Sentinel-2和一些自己的GPX文件等）对比，选取了一些长期不变的点\n",
    "#path_points = '/mnt/2thdd-btrfs/rs/USGS_Declass/DS1015-1021DF073/DS1015-1021DF073_c.tif.points'\n",
    "path_points = '/mnt/2thdd-btrfs/rs/USGS_Declass/D3C1215-401419A011/D3C1215-401419A011_e.tif.points'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "0755fc7e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "30.08455, 120.56050, 5, 16433, -17237\n",
      "30.08066, 120.55324, 6, 17228, -16622\n",
      "29.99803, 120.56361, 8, 14563, -5574\n",
      "30.05613, 120.45384, 9, 28652, -12034\n",
      "30.05347, 120.45612, 8, 28340, -11704\n",
      "30.03196, 120.45630, 15, 27982, -8784\n",
      "30.00135, 120.53514, 6, 18020, -5647\n",
      "30.07184, 120.41992, 8, 32957, -13729\n",
      "30.12122, 120.60107, 7, 12294, -22628\n",
      "30.09488, 120.43264, 10, 31790, -17001\n",
      "30.05497, 120.66658, 6, 3449, -14616\n",
      "30.10485, 120.39856, 7, 36010, -17904\n",
      "29.98058, 120.65131, 6, 3855, -4364\n",
      "30.00076, 120.43992, 38, 29444, -4307\n",
      "30.04988, 120.46144, 8, 27644, -11285\n",
      "30.12739, 120.40110, 12, 36037, -20964\n",
      "30.10306, 120.50965, 5, 22764, -19073\n",
      "30.07040, 120.47227, 12, 26686, -14204\n",
      "30.02731, 120.53981, 4, 17908, -9235\n",
      "29.99726, 120.42892, 49, 30722, -3682\n",
      "30.01593, 120.42640, 51, 31305, -6205\n",
      "29.98720, 120.41517, 128, 32218, -2124\n",
      "29.98738, 120.42289, 193, 31290, -2270\n",
      "29.98808, 120.38073, 151, 36405, -1800\n",
      "30.03474, 120.45713, 10, 27918, -9164\n",
      "30.02654, 120.43360, 11, 30614, -7748\n",
      "30.02660, 120.39368, 177, 35413, -7256\n",
      "29.98942, 120.45608, 81, 27321, -2965\n",
      "30.00947, 120.45178, 63, 28157, -5658\n",
      "30.08213, 120.45197, 7, 29289, -15527\n",
      "30.06716, 120.39356, 56, 36044, -12756\n",
      "30.06001, 120.63088, 3, 7708, -14833\n",
      "30.01763, 120.62717, 7, 7378, -9081\n",
      "30.13555, 120.44241, 10, 31254, -22565\n",
      "30.09650, 120.66382, 5, 4527, -20110\n",
      "30.09465, 120.66554, 8, 4294, -19889\n",
      "30.11318, 120.67659, 4, 3354, -22470\n",
      "30.10925, 120.68760, 12, 2000, -22076\n",
      "30.10772, 120.69552, 16, 1059, -21977\n",
      "30.09555, 120.61877, 9, 9772, -19433\n",
      "30.03746, 120.51172, 10, 21430, -10255\n",
      "30.03270, 120.50853, 10, 21736, -9568\n",
      "30.05593, 120.58114, 7, 13495, -13653\n",
      "30.05934, 120.69211, 5, 561, -15519\n",
      "30.01531, 120.64773, 10, 4920, -9039\n",
      "30.00219, 120.40846, 225, 33261, -4111\n",
      "30.01040, 120.40906, 202, 33315, -5245\n",
      "30.07165, 120.52082, 7, 20920, -15000\n",
      "30.11197, 120.58314, 7, 14243, -21170\n",
      "30.11100, 120.58091, 11, 14493, -21018\n",
      "30.10705, 120.58323, 8, 14143, -20529\n",
      "30.07521, 120.58737, 8, 13098, -16329\n",
      "30.08566, 120.54288, 9, 18530, -17162\n",
      "30.08474, 120.54726, 11, 18000, -17097\n",
      "30.06353, 120.46619, 7, 27289, -13197\n",
      "30.06401, 120.44216, 8, 30175, -12957\n",
      "30.06091, 120.44518, 10, 29764, -12572\n",
      "30.07572, 120.41728, 8, 33333, -14222\n",
      "30.03908, 120.49031, 8, 24009, -10197\n",
      "29.99101, 120.47639, 5, 24893, -3446\n",
      "29.98795, 120.47208, 7, 25362, -2966\n",
      "29.97991, 120.48364, 1, 23842, -2012\n",
      "29.98088, 120.66274, 116, 2517, -4563\n",
      "29.97532, 120.58484, 8, 11631, -2754\n",
      "29.98376, 120.57815, 11, 12575, -3817\n",
      "29.98120, 120.57301, 9, 13140, -3399\n",
      "30.01409, 120.48345, 9, 24428, -6702\n",
      "[    30.04795846    120.5175076      29.40102109  20973.17184434\n",
      " -11716.46344994]\n",
      "67\n"
     ]
    }
   ],
   "source": [
    "osra, res = read_points(path_points)\n",
    "osrgeo = osr.SpatialReference()\n",
    "osrgeo.SetWellKnownGeogCS(\"WGS84\")\n",
    "ct_geo2prj = osr.CoordinateTransformation(osrgeo, osra)\n",
    "ct_prj2geo = osr.CoordinateTransformation(osra, osrgeo)\n",
    "data = []\n",
    "for mapX,mapY,sourceX,sourceY,enable,dX,dY,residual in res:\n",
    "    lat, lon = ct_prj2geo.TransformPoint(mapX, mapY)[:2]\n",
    "    h = dem.get_geopoint_cubic(lon, lat)\n",
    "    data.append((lat, lon, h, sourceX, sourceY))\n",
    "    print(f'{lat:.5f}, {lon:.5f}, {h:.0f}, {sourceX:.0f}, {sourceY:.0f}')\n",
    "data = np.array(data)\n",
    "print(data.mean(axis=0))\n",
    "print(len(data))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 246,
   "id": "f40eac7d",
   "metadata": {},
   "outputs": [],
   "source": [
    "pl=(data[:,:2]-(30.05, 120.5))\n",
    "plh=np.concatenate([pl, data[:,2:3]], axis=1)\n",
    "#plh=pl\n",
    "xy = data[:,-2:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 247,
   "id": "5133706e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.00000000e+00  1.65570159e+04 -1.19354123e+05 -1.81760158e+03\n",
      "   1.35945588e+04  5.17128457e+03 -3.08640316e+03 -1.87303646e+03\n",
      "  -1.26359432e+04  1.18564360e+04 -2.68082490e-02]\n",
      " [ 0.00000000e+00 -1.35503229e+05 -1.30097399e+04  1.47300284e+04\n",
      "   6.60923803e+03  3.94938322e+01  1.05064602e+04  1.09189307e+04\n",
      "   1.26825453e+04  4.49654544e+03 -1.32634780e-01]]\n"
     ]
    }
   ],
   "source": [
    "pfX = PolynomialFeatures(degree=3)\n",
    "pfX.fit(pl)\n",
    "def mytransf(x):\n",
    "    #x = pfX.transform(x)\n",
    "    # remove z^3\n",
    "    #return np.concatenate([x[:, :19]], axis=-1)\n",
    "    # remove z^2, xz^2, yz^2, z^3\n",
    "    #return np.concatenate([x[:, :9], x[:, 10:15], x[:, 16:18]], axis=-1)\n",
    "    # remove z^2, x^2z, y^2z, xz^2, yz^2, xyz, z^3\n",
    "    #return np.concatenate([x[:, :9], x[:, 10:12], x[:, 13:14], x[:, 16:17]], axis=-1)\n",
    "    xpl = pfX.transform(x[:,:2])\n",
    "    return np.concatenate([xpl, x[:,2:3]], axis=-1)\n",
    "X = mytransf(plh)\n",
    "model=linear_model.LinearRegression()\n",
    "model.fit(X, xy)\n",
    "print(model.coef_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 248,
   "id": "05d36b3d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1.58675689 1.85915527 3.01901287 1.61363511 1.87755909 7.98008725\n",
      " 2.51729369 3.16338206 5.2953873  2.69439071 3.68775675 5.52887314\n",
      " 2.13969136 3.77608942 0.94765847 0.63566668 2.58847367 6.48129373\n",
      " 9.01092092 3.85291144 9.70656705 5.26233175 7.87100413 1.27111921\n",
      " 6.9165604  1.92352466 3.15258421 6.18156146 0.85662591 0.21124665\n",
      " 6.04541544 6.98728056 2.99950974 3.37017933 1.60181305 1.4751076\n",
      " 3.14250124 6.6877216  3.71466735 3.09275338 1.55372744 5.04015909\n",
      " 2.66155811 2.86993487 4.41458552 1.7573647  2.7071351  2.74472254\n",
      " 7.95713121 9.10987167 4.4719062  3.52490432 7.62245475 4.47875505\n",
      " 5.54096932 5.11457583 1.6663186  4.28594511 5.43085665 2.32467625\n",
      " 3.50920426 6.29846569 8.24148642 3.97182736 2.94223422 6.09498456\n",
      " 3.1626378 ]\n",
      "9.706567047510687 3.612774126154679\n"
     ]
    }
   ],
   "source": [
    "err=model.predict(mytransf(plh))-xy\n",
    "errx=err[:,0]\n",
    "erry=err[:,1]\n",
    "erra=np.sqrt(errx**2+erry**2)\n",
    "rms=np.sqrt(np.sum(err**2)/(np.size(err)-np.size(model.coef_)))\n",
    "print(erra)\n",
    "print(np.max(erra), rms)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 252,
   "id": "3717de54",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAD4CAYAAAAHHSreAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/SrBM8AAAACXBIWXMAAAsTAAALEwEAmpwYAAA3mklEQVR4nO3deXxU5bnA8d8zS1aWsAQJmyCLCAioAcG9UBWXqrXUamur1ta26tVWu1l7bW1rq95bt7baS6UttlbErVJcKioWUQGDyL6FTQKBJED2Zbbn/jEnGDITMkkmmUzyfD+ffD6ZM+c95zmZJM95l/O+oqoYY4wxDbkSHYAxxpjOx5KDMcaYCJYcjDHGRLDkYIwxJoIlB2OMMRE8iQ6gJfr376/Dhw9PdBjGGJNUVq1aVaKq2S0pk1TJYfjw4eTl5SU6DGOMSSoisrulZZptVhKRNBFZKSJrRGSDiNzrbB8hIitEJF9EnhWRlChl+4nIEhGpFJHfN3H8hSKyvqWBG2OMaT+x9DnUATNUdRIwGZglItOAB4CHVXUUcBi4MUrZWuC/ge9HO7CIXAlUtiJuY4wx7ajZ5KBh9f/Avc6XAjOA553t84AropStUtVlhJPEUUSkB3AH8KtWRW6MMabdxDRaSUTcIvIxUAQsBrYDpaoacHYpAAa38Ny/BH4LVDdz7ptEJE9E8oqLi1t4CmOMMa0RU3JQ1aCqTgaGAFOBsW05qYhMBkaq6ksxnHuOquaqam52dos6240xxrRSi55zUNVSYAkwHcgSkfrRTkOAvS041HQgV0R2AcuAMSLyTktiMcYY035iGa2ULSJZzvfpwPnAJsJJYraz23XAy7GeVFWfUNVBqjocOAvYqqrntShyEyEUUu56cS3jf/Y6Fz26lD2HjtliZ4wxTYql5pADLBGRtcCHwGJVXQT8CLhDRPKBfsBcABG5TER+UV/YqR08BFwvIgUiMi7O12AcL6/Zyz9X76OqLsiW/RXcPn91okMyxiSpZh+CU9W1wClRtu8g3P/QePtCYGGD18ObOf4uYELzoZrmFJbV4guGAAhp+LUxxrSGza3UhVw8IYd0r5t0r4t0r5uvnzki0SEZY5JUUk2fYY5teP9MXv/u2SzdWsLw/hmcMbJ/okMyxiQpSw5dzJA+GXz59GGJDsMYk+SsWckYY0wESw7GGGMiWHIwxhgTwZKDMcaYCJYcjDHGRLDk0EKqylMf7GJ7sS1DYYzpuiw5tFCtP8Qji7fx7/X7Ex2KMca0G3vOoYXSU9ws+cF59Ey1H50xpuuy/3Ct0Dvdm+gQjDGmXVmzkjHGmAiWHIwxLaKqbD1Qwfq9ZQRDmuhwTDuxZiVjTMwOlNfy1bkr2HOoBpdAZqqHv94wlXGDeiU6NBNnsawElyYiK0VkjYhsEJF7ne0jRGSFiOSLyLMikhKlbD8RWSIilSLy+wbbM0TkFRHZ7Bzz/vheljGmPdw+fzXbiyqp8Qep8gUpqqjj+r+sJGQ1iC4nlmalOmCGqk4CJgOzRGQa8ADwsKqOAg4DN0YpWwv8N/D9KO/9r6qOJbyQ0JkiclEr4jfGdJAaX5APdx0m2CgPVNUF2LS/PDFBmXbTbHLQsPonvrzOlwIzgOed7fOAK6KUrVLVZYSTRMPt1aq6xPneB3wEDGnlNRhjOoDLBS6J3B5SSPW4Oz4g065i6pAWEbeIfAwUAYuB7UCpqgacXQqAwa0JQESygM8BbzXx/k0ikiciecXFxa05hTEmDlI9bi4+OYdUz6f/NtwuGN4/g1EDeiQwMtMeYkoOqhpU1cmE7+6nAmPjcXIR8QDPAI85a1JHO/ccVc1V1dzs7Ox4nNYY00r3XzmRiybk4HULbpcw7YR+zPt6xFLypgto0WglVS0VkSXAdCBLRDxO7WEIsLcV558DbFPVR1pR1hjTwdJT3Dxy9WQenD2RkCppXmtO6qpiGa2U7TT9ICLpwPnAJmAJMNvZ7Trg5ZacWER+BfQGvtuScsaYxEvxuCwxdHGx1BxygHki4iacTBao6iIR2QjMd/7JrwbmAojIZUCuqt7jvN4F9AJSROQK4AKgHLgb2Ax8JCIAv1fVJ+N4bcYYY1qp2eSgqmsJDzdtvH0H4f6HxtsXAgsbvB7exKGjjHswxrQnXyDEgfJaMlLc9OuRmuhwTCdmT0gb08WV1/p5efU+XlpdwPq9ZbhcQigEY3N68vPLxnPqsD6JDrFL2lVSRZUvwLicXjitI0nFkoMxXZSq8vg723nsrW34gyGOPMTsPMW2tqCML/9pOX+78XSmDO+buEC7oD/+ZzuPvLkVQTh3TDZPXHtq0iUIm3jPmC5qzrs7+P3b+dQFGiSGRmr9Ie56cV3HBtYNPPrmNmr9IWr8Qd7ZUkTB4ZpEh9RilhyM6YKqfQEeWbyNGn+w2X0LDlezs6SqA6LqPnpnfLrmiwI905pupPEHQ+wqqeJgZV0HRBY7a1YypgtavPEArhhv/bxuF6XVPiCzXWPqTp78Wi63zV9NVV2Aey4dR1ZGxLykAPzlvZ08vHgrgZASCCqnDMvi4S9NZlBWegdHHMmSgzFdUHFFHb5AKKZ9fYEQx/ezxBBPEwb35u07zzvmPn9fvpsHX99yVO0ub/dhPv/4e/znB59J+HMk1qxkTBfUr0cKKe7m/7wFOO/EbPpmRr+zNe0jFFIeXrw1otkvGFIqawMsWluYoMg+ZcnBmC7o/HEDCWrzaywM7J3G/VdO7ICITEOHq31U1AaivlflC5K36xAQHoZ854KPuejRpTz4+mYCwdhqg/FgzUrGdEE9Uj3cdM4J/GnpziY7pdM8Ll67/ewm28NN+8lMbfpfr9ctHNcrDYBbn/6I5TsO4gsqu0qqcLuEOy84sUNitJpDjMpq/Ow+WGVr5pqk8b3PjuGr048n1eM6aprthl7+eF8HR2UA0rxuLpowEK878tkHlwizTwsvb5O3+zA+57mUGn+Id7Z03LIFVnNohi8Q4kcvrOWVdYW4RUhPcfPbL07iM2MHJDo0Y45JRPjJxSfxjbNH8NJHe3lo8VbqGnRS1wZC7D5oQ1gT5RdXTGDLgQr2HKqmyhc8ksB/c+XJDO2bAcAJ2Zls2ldBUJUUj4uTB/fusPhEY2iX7Cxyc3M1Ly+vQ8/5y0UbeXrFbmr9n/5RpXldvHnHuQzpk9GhsRjTFtf/eSXL8ksIOLXfdK+bB2ZP5LJJgxIcWfcVCilLtxWTt+swWRleLps8iAE90468v7e0hu/8fRU7i6uYNrIfj3xp8jGbpJoiIqtUNbdFZSw5HNuEn/2byrqjO468buH2maO5dcboDo3FmLYoqazj2idXsLMk3Dx67bTj+dnnxiXdtA6m5VqTHKxZqRnRxooHQ0q1r/knT01iHCiv5V9r9jEupxdnjOqf6HA6jf49Unnt9rMprqwj3eumZ5q3+UKm27Lk0IwZY7N5c1PRkao4hBc6uWhCTgKjMk0JBEN87nfLKK324XIJc6+bwpmWII4QkaOaLYxpio1WasYvrziZoX3TyUx1k5niJtXj4tvnjOTkIR3XMWRiV1kX4FCVD19QCYaUDfvKEh2SMUmp2ZqDiKQBS4FUZ//nVfVnIjICmA/0A1YBX1VVX6Oy/YDngSnAX1X11gbvnQb8FUgHXgVu107YAZLdM5W37jiP5TsOUlRRx5QRfRncCeY9MdFlZaRw6cRBLFq7j15pXi6daJ2txrRGsx3SEu6tylTVShHxAsuA24E7gBdVdb6I/BFYo6pPNCqbSXgVuQnAhEbJYSVwG7CCcHJ4TFVfO1YsieiQNsnpUJWPnmkevDFMIWFMV9eaDulm/3I0rNJ56XW+FJhBuFYAMA+4IkrZKlVdBtQ2CjQH6KWqy53awlPRyhvTWn0zUywxGNMGMf31iIhbRD4GioDFwHagVFXrx3gWAINbcN7BTpl6TZYXkZtEJE9E8oqLO+7pQGOM6c5iSg6qGlTVycAQYCowtj2DanTuOaqaq6q52dnZHXVaY4zp1lpU71bVUmAJMB3IEpH6Du0hwN4WHGqvU6ZeS8sbY4xpR80mBxHJFpEs5/t04HxgE+EkMdvZ7Trg5VhPqqqFQLmITHM6vL/WkvLGmOSmqry2rpB57++ipJMtjxnNxn3lnP3A29y54ONEh9JhYqk55ABLRGQt8CGwWFUXAT8C7hCRfMLDWecCiMhlIvKL+sIisgt4CLheRApEZJzz1s3Ak0A+4T6MY45UMvBefgn/WPEJZdX+RIdiTJvc/9pm7nxuDb9+dRMXPfouVXXR1zboLF5dV0jB4Rpe+GgvtTGsy90V2NxKSWLhx3v50QvrUJSBvdJY8v3zbE4ck7TOvP9t9pbWAOG1J566cSqnDuuT4Kiatq+0hp++tJ4TsjNI9bj52hnDj6y5kAzaZSir6RyW5ZdQ4w9S6w+x53BNxGSAxiST047vc9QaEyM6+RrWg7LS+fMNU1i6rYTH39nO/a9tTnRI7c7mVkoSXzh1CP9aU4gI5B7fhx6tmLbXmM7iwdkTGd4vg8KyWm44cwR9kmQN68smDeJP7+5g1oSBiQ6l3VmzUhIpLKthf1ktE4dk4XZZk5IxJjY2ZXcXl9M7nZzeNq+TMab9WZ+DMcaYCJYcjDHGRLDkYIwxJoL1ORgTR6GQ8p9txXy06zDpqW4umpDDiP6de5imMdFYcjAmTvaW1nDNnOUcrKyjyhfE4xIefXMbl08exP1XTsRlI8xMErFmJWPiQFW59snl7D1cQ5UvPL1CIKTUBUL8a00hT/xne4IjNKZlLDkYEwcfbD9IUXkdwSjPDdX4g/z2jS1c+PBS5q/8JAHRGdNylhyMiYOPC0qPOSFbSGHLgQru/ddGnv0wtgRxuMpHja97TPJmOh9LDsbEQWaKB08My5LW+IM89cHuZvdbtfswp//6Lc568G2qfTaPlul4lhyMiYOWzLWT4mn+z67gcDUuF1TUBKiqs9qD6Xg2WsmYODiuVxo3nDmcee/vpuYYzUtpXhe3zxzd7PE+N3EQwZAyKCud7J6p8QzVmJhYcuimthdXkupxMaRPRqJD6TJ+NGssQ/tm8Nhb2zhU5UM1PDX1l08fxtKtxfiDyjVThzF9ZL9mj+VyCVeeOqTZ/YxpL80mBxFJA5YCqc7+z6vqz0RkBDCf8Cpwq4CvqqovSvm7gBuBIHCbqv7b2f494BuAAuuAG1S1Ni5XZZpVUlFHzzQvdN71VZKOiPCV04/ny1OHUVrtJ9XrIiMl/Cd2+eTBEfv7AiGeX1XA0L7pnD06u6PDNeaYYulzqANmqOokYDIwS0SmAQ8AD6vqKOAw4QRwFGdJ0KuB8cAs4HERcYvIYOA2IFdVJwBuZz/TQU4/oR/jBvVKdBhdkojQJzPlSGJoyh+W5HPvvzbwzXl5bDtQ0UHRGRObZpODhlU6L73OlwIzgOed7fOAK6IUvxyYr6p1qrqT8HrRU533PEC6iHiADGBfay/CmGSU5nXhEkDAG8NIJ2M6Ukx9DiLiJtx0NAr4A7AdKFXV+jF2BUBkvTm8bXmD1wXAYFX9QET+F/gEqAHeUNU3mjj3TcBNAMOGDYslXGOSwk3njGRE/0wGZaUz3OZfMp1MTLcrqhpU1cnAEMJ3/mPbclIR6UO4VjECGARkisi1TZx7jqrmqmpudra1y5quw+0SZk3IYeKQrESHYkyEFtVlVbUUWAJMB7KcJiEIJ429UYrsBYY2eF2/32eBnaparKp+4EXgjJaFnnzqAkHmvruDu19ax7/W7COZlmg1xnQvsYxWygb8qloqIunA+YQ7o5cAswmPWLoOeDlK8YXAP0TkIcI1hNHASiAETBORDMLNSjOBLr04dCikXPvkCtYWlFEXCPHiR3vZsr+C7194YqJDM93YtgMVzP9wD30yvNxw5ggyU210uwmL5TchB5jn9Du4gAWqukhENgLzReRXwGpgLoCIXEZ4FNI9qrpBRBYAG4EAcIuqBoEVIvI88JGzfTUwJ94X15ls2l/Ohn3l1AVCQHgahf9bup07LxiDiE3lbDrenkPVXPGH96jyBUnxuHhjwwFevvVM+300QAzJQVXXAqdE2b6DT0ceNdy+kHCNof71fcB9Ufb7GfCzFsabtIIhpfGfnGr4y/4WTSK8l19CyGna9AVCbNpfTlmNn6yMlARHZjoDGz/XQU7K6cWgrHRS3OFMkO51M/u0IbYAjEmYYX0zoMEtS6rHTQ9rVjIOSw4dxOt28cLNZ/CV04/nnDHZfPezo7l22vFsL65svnAnVxcIsre0xqaXTjJnjOrPzeeNpGeah0FZafzlhikxzSxrugdJphEzubm5mpfXNfqtb5+/mjc2HEBRvnn2Cdx5QfJ1TIdCykOLt/Dn93Y5TWTKF04bwj2fG0eqx53o8Jqlqjy8eCuvrivk0omDuP2zo6293XRJIrJKVXNbUsZuExKgqLyW19btp8YfpNYf4vF3tiflsNb7X9/M3GW7qPYFw9cSCPHCqgLuWLAm0aHFZFl+CU8u20l+cRVz3t3B+9sPJjokYzoNSw4JkJ7iPqoTukeqJ+nuWKt9AZ76YFfE9NS1gRCLNx6gsKwmQZHFrnEzmDWLGfMpSw4J0DPNy6NXTya7ZyqDs9J58roW1fY6hU8OVeNxRf/1SfW42Hqg8/elzBg7gBljB5CZ4mbmSQP4zNgBiQ7JmE7DhiY0snjjAR56YwsFh2sYNaAHP5w1Nqb591tq1oQcZk3IiftxO8qAnmn4gqGo7wWCSk7vtA6OqOU8bhe///KpiQ7DmE7Jag4NvPRRAf/1zEds2l9BRV2A1XtKueGvK1m2rSTRoXU6fTNTOHdMNl730c1hboGR2ZmMOa5ngiIzxsSDJQeHqvLrVzdT6z/6brjWH+LXr25KUFSd22+vmsQpw/qQ5nXRI9VNRoqbUQN68ufrpyQ6NGNMG1mzkqO8NkBpTcRCdgBsK7KFWKLpleZlwbems2V/BduKKhjaJ4OJQ3onXee6MSaSJQdHZoobr9uFPxg5YiW7R3Is8F7rD7JwzT42F5Zz6vF9mDV+YFwfavIHQ+w+WEX/HqlHTbFw4sCenDjQmpGM6UosOTg8bhc3nDmcPy87enhmutfNrTNGJzCy2FT7Anzud8vYV1ZLjS/I/A/38NT7u3nmpmm44zBFx6K1+7jrxXWEQoo/pFx6cg4PzJ5oK5gZ00XZX3YDd5x/ItdOG0aa10W6101mqofbZo7imqlDmy+cYM/lFbCvtPbIWP1qX5D1+8pYsrmozcfedqCC7z+3horaAFW+IL5AiFfXF/LI4q1tPrYxpnOymkMDbpdw9yXjuPOCEzlY5SO7RyopnuTIn6s/ORzxQFqNL8imwnI+O+64Nh37uVUF+AORHfX/WPkJP5jVpkUBjTGdlCWHKNK8bgZnpSc6jBaZPDSL1zfsP2q0VXqKm5NyerX52DW+IKEos3v4g8k35YdJPqrK+9sP8s6WIlI8Li45eRDjBrX999ocW7O3xSKSJiIrRWSNiGwQkXud7SNEZIWI5IvIsyISdRJ4EbnL2WeLiFzYYHuWiDwvIptFZJOITI/fZXU/V00ZyuCsdNK94QnvMlLcTBjUm/NObPu625dOzCHNe/REel6XcMH4ttVIjGlOjS/IF554n28+lcef3t3JE+9s58on3uMHz61JyvnIkkksbSZ1wAxVnQRMBmaJyDTCS4U+rKqjgMPAjY0Lisg44GpgPDALeNxZUQ7gUeB1VR0LTALsYYI2yEjx8MptZ3Pv5eP5+pnDeXD2RPr3TGHMT1/jubw9bTr26Sf047ozhpPqcZGZ4iYzxc0J2T2459JxcYremOjue3UTG/aVU+30pYU03KS5aG0hL3wUbdl6Ey+xrASnQP1EOV7nS4EZwJed7fOAnwNPNCp+OTBfVeuAnSKSD0x1lhg9B7jeOYcPiP6QgYlZmtfNVbmfdp7/8Pm1hBReX7+fL+a2rVP9xxeN5SunD+OjTw6T0zudKcP72PMMpl35gyGeX7XnyNK6DdX4g8xZup3Zpw1JQGTtq6zGz4Ovb6a02s+dF4zhhOweCYkjpj4H525/FTAK+AOwHShV1YCzSwEwOErRwcDyBq/r96sBioG/iMgk59i3q2pVlHPfBNwEMGzYsFjCNY5Hr57MwjX7uOP8+KwVMbRvBkP7ZsTlWMY0p6ouQDBaZ5fjQHldB0bTce549mOWbi0moMrKXYdY+ZOZCbkRi2kojqoGVXUyMITwutFtHaLiAU4FnlDVU4Aq4MdNnHuOquaqam52dtvbz7uT88cN5HfXnMqI/pmJDsWYFuuZ5o3o62poeL+ueaOyo6QKf0hRhYOVdcdMkO2pReM0VbUUWAJMB7JEpL7mMQSI1gC4F2jYnlG/XwFQoKornO3PE04WxhgDhIeWf/3MEaR7I/9NpXvd/FcSPJzaGt+dOZpUj4tUj4uvTR+esKVbm21WEpFswK+qpSKSDpxPuDN6CTAbmA9cB7wcpfhC4B8i8hAwCBgNrFTVoIjsEZETVXULMBPYGJcrMi1S7QuQ5nHjisNT1KbzqgsEeWdLMQN7pTFpaFaiw4nZbTNH88mhal5dV4gIuEQIhpRbZ4xq8/M7ndXlpwxm+qh+1PpCDEtg7SiWPoccYJ7T7+ACFqjqIqdTeb6I/ApYDcwFEJHLgFxVvUdVN4jIAsL/+APALapa/6TWfwFPO0NgdwA3xPXKTLMWbzzAt/+2irE5PVn0X2dZB3MXdv2fP2RtQSkhhQe+cDKXTY7WRdj5uF3Cw1+azHc/O5pl+SWkuF3MGDuAfgma78wXCPG7t7exrqCM88Zmc9304e3ydzOgZ+LXQ4lltNJa4JQo23cQ7n9ovH0h4RpD/ev7gPui7PcxkHxLoHUhqz85jAhs2V9BMKR43JYcuiJVZfnOg9Q/FvD6hv1JkxzqHd8vk+P7Jb7v7IfPrznysOmKnYeo84f41rkjEx1Wu7AnpLuxmz8ziowUN6ce3ydh7Zqm/YkIM8YO4IPtBwmp8vlTut7wz47y9uaiI7MQ1PiDvLKu0JKD6Xp6pHq4dcZodpZUMf03b5GR4mbBt6YnrMpu2s//XXsaK3cdYkDPVEYNsOnVW2tI3ww2FZajCl63MCpBzyB0BLtdNLy2vpCi8lr2Hq5hWb4tidoVedwuzhjZ3xJDGz3xlfDQcLdLmDQki599bnyiQ2o3VnMwXDwhh6eXf0JmqpuzR9uzJMY05fh+mbx953mJDqNDWHIwDO+fyXs/npHoMIwxnYg1K5kmhUJKIBg5r40xpuuz5GCadNX/fUDur97EbwnCmG7HmpVMk4b0SccfDOGyh+OM6XYsOZgmPXJ1xLOPxphuwpqVjDHGRLDkYIwxJoIlB2OMMREsORhjjIlgycEYY0wESw7GGGMiWHIwxhgTodnkICJpIrJSRNaIyAYRudfZPkJEVohIvog866zoFq38Xc4+W0TkwkbvuUVktYgsis/lGGOMiYdYHoKrA2aoaqWIeIFlIvIacAfwsKrOF5E/AjcCTzQsKCLjgKuB8YTXkH5TRMY0WCr0dmAT0Cs+l2OSUbUvwH2vbOLlj/cBcMXkQdx9yTjSU9wJjsyY7qvZmoOGVTovvc6XAjOA553t84ArohS/HJivqnWquhPIx1laVESGAJcAT7blAkzyu/npj3huVQGVdQEq6wI8t6qAm59eleiwjOnWYupzcJp/PgaKgMXAdqBUVQPOLgVAtEVpBwN7GrxuuN8jwA+BY87qJiI3iUieiOQVFxfHEq5JIvtKa/hg+0F8gU9/DeoCId7ffpDCspoERmZM9xZTclDVoKpOBoYQvvMf25aTisilQJGqNnt7qKpzVDVXVXOzs20hmq7mUJUPb5T1qz1u4WClr83HD4W0zccwpjtq0WglVS0FlgDTgSwRqe+zGALsjVJkLzC0wev6/c4ELhORXcB8YIaI/L1FkZu4OlTl4/X1hby58QDVvkDzBeJkzHHRl610iTT5Xqze2nSAUXe/ytf/8mGbjmNMd9Rsh7SIZAN+VS0VkXTgfOABwkliNuF/7tcBL0cpvhD4h4g8RLhDejSwUlU/AO5yjn8e8H1VvbbNV2NaTFV58PUtzH1vJynOHXwwpPz8svF8acrQZkq3XYrHxSNfmsytz3x01PZHvjSZFE/bRlqv3HkIEWHlrkNtOo4x3VEso5VygHki4iZc01igqotEZCMwX0R+BawG5gKIyGVArqreo6obRGQBsBEIALc0GKlkOoG/L9/NX9/fhS8QOqrd/+cLNzCsbwbTR/Zr9xg+O+44lv7gM/x7w34ALpwwkAE909p83FtmjKJnmqdDrsGYrkZUk6dNNjc3V/Py8hIdRpcy5b43Ka6oi/re9JH9eOab0zo4ImNMvInIKlXNbUkZW+ynG6us8zeZGAC2HajowGjabs+han6+cAOrdh9mWL8MfnrJOKaO6JvosIxJSjZ9Rjd278KNx3w/u0dqB0XSdrX+IJ9//D2WbCmitMbP2oIyrvvzSvKLkivBGdNZWHLopkoq61i4Zt8x9zlnTPIMHX57cxE1/iANR67WBYL8bfkniQvKmCRmyaGben/7QTxuOeY+lbX+Doqm7arqAjTuPgtpcl2DMZ2JJYduSlXDk6AcQ4hjJ4/OZMbYAYQaZYd0r5vPnzKkQ85fUetnbUEpG/eVEwge86F/Y5KCdUh3U1NH9MV/jKeHM1PdnHdi8jQr9euRyhPXnsYdz35MrT+Eotw2cxRnje7fruc9XOXjl4s28sq6QlLcLkIoXpeLb517At86ZyQuV/IkWGMasuTQTeX0Tmfm2AG8vbmIusDRd7ougaz0FGaOHZCg6FrnMycO4MO7P0thWS39e6S2+6yupdU+Lv3dMooqavEHtcHPMchjb+WzcV85j11zCiKWIEzysWalbuyhqyZz+oi+pHld1N/gZqa4Gdo3gwXfno4nypxHnZ3H7WJo34wOme77t29sPZIYGqvxB3lrcxHL8kvaPQ5j2oPVHKLYUVzJip2HOHlwbyYM7p3ocNpNeoqbp248nY37ynlz0wH8gRBTRvTlrFH9rTmkGXWBIC98VBA1MdSr9gWZs3QHZ49OnuY5Y+pZcmhk8/5yrnz8fVQVBeZ8NbfTDOkMhpRf/GsDO0uq+M0XJjI4Kz0uxx03qBfjBtl6Sy1RVN70w4MNbdlvz1mY5JR87QbtbPGGA9T6g9T4Q9T6Qzy/qiDRIR2xYsdBFuQVsCy/hN+9tS3R4XRrqV4XgWPUGuqlee1PzCQn+81tZPRxPUn1hNur071uJgzuPHfUJ2T3wOMWPG4Xp9u0EAk1oGcaQ/seu+aW4nHxuUmDOigiY+LLmpUauXD8cfxo1om8sq6QaSf048azTkh0SEcM7J3Gsh/NoLIuELcmJdN6d15wIncuWEONP/pEwx6X8LXpwzs2KGPixJJDIyLC9WeO4PozRyQ6lKh6p3vpne5NdBgGuPjkHPKLKnl8ST7+UIj6Z9/SvC7cIsy9fgrH9Wr71OPGJIIlB2Pa4LaZo7lowkD+/N4uVu0+RIrbxaUTc7hqyjD6ZqYkOjxjWs2SgzFtNPq4nvzmypMTHYYxcdVsh7SIpInIShFZIyIbROReZ/sIEVkhIvki8qyIRL1NEpG7nH22iMiFzrahIrJERDY6x7w9vpdljDGmLWIZrVQHzFDVScBkYJaITCO8jvTDqjoKOAzc2LigiIwDrgbGA7OAx53lRgPAnao6DpgG3OLsa4wxphNotllJw+uIVjovvc6XAjOALzvb5wE/B55oVPxyYL6q1gE7RSQfmKqqHwCFzvErRGQTMJjwWtPd2po9pfxnazG1/iCjj+vBRRNySPO2/1QQxhjTUEx9Ds7d/ipgFPAHYDtQqqoBZ5cCwv/cGxsMLG/wOmI/ERkOnAKsaOLcNwE3AQwbNiyWcJPSJwer+cZTH7LnUA11gfCiNZkpbu5+aT0/veQkvnz68YkO0ZgOtbe0hj2Hqpk0JKtD5soyR4spOahqEJgsIlnAS8DYeJxcRHoALwDfVdXyJs49B5gDkJub2/wjqUnoQHktl/9hGWU1/qNWMqvyhcfP/3LRJlThK9MsQZju4f38Em6cl4fLFV6u9tXbzyYjxcbPdKQWPSGtqqXAEmA6kCUi9Z/WEGBvlCJ7gaENXh/ZT0S8hBPD06r6YsvC7loeeXMbFbUBmlpeocYf5L5XN1Hji/6wlTFdzdz3dlLjD1JVF6S4so6VOw+1+ZhVdQGKKmrDC12ZZsUyWinbqTEgIunA+cAmwklitrPbdcDLUYovBK4WkVQRGQGMBlZKeIL7ucAmVX2ozVeRxGr9QV5aXUDgGAvv1HtlXWEHRGRM4o0Z0OPIvFShkDK0b0abjvfbN7Yw6d43OPuBJVzw8FKKK2KbOLE7i6WelgPMc/odXMACVV0kIhuB+SLyK2A14X/2iMhlQK6q3qOqG0RkAeGO5gBwi6oGReQs4KvAOhH52DnPT1T11bheXRIoKq/DFcNiMNW+IJv3R215M6bL+e75Y/AHlQ37yrnhzOGMzO7R6mMt2VLEk+/uJBBSAiFlZ0kV33v2Y/7+jdPjGHHXE8topbWEO4wbb98BTI2yfSHhGkP96/uA+xrtswySaIHiduT1CMEYag0CpHmSu1Nuf1kt185dwa6SKi4cP5DHrjkFt60bYaJI9bj56aXxGd2+YW8ZdYFPm2QDIWX9vrK4HLsrs1lZE2xgrzT69Wh+moX0lORa0zmaX76ykR3FlQRCypItRdZMZjrEsH6ZRw0HF4FhbWym6g4sOSSYiPCdc0eSfoxnGURgQM9UTju+TwdGFn9VDTrdVZVa62A3HeDSk3M4Z0w26V43PVM99ElP4aGrJic6rE7PxoZ1Al85/XiWbi1hWX5JxPTPLoEeqR7mXj8l6Req/8GsE1m1+zCBkDK0bzqXTMxJdEimG3C5hCe+ciqbCiuoqPUzblAveqbZzMbNkWQa1pWbm6t5eXmJDuOYquoC7DpYxagBPY4sGhSLYEj52we7+OPSHZRW+3CLEAgpl07M4Xvnj2FIn65RDa72BSiuqGNwVjoet1VcjekIIrJKVXNbVMaSQ/wUldcy69F3qfUHGdCzdQ/uqCoFh2vwBUPk9E6zB3+MMW3WmuRgt25xtGRLETW+INW+ICWVdXz8SWmLjyEiDO2bwcjsHu2aGA5W1hGoX53GGGMaseQQR+MH9UZR3C5BNbzmc2e0r7SG0371Jvcs3JDoUIwxnZS1WcTRhMG9efobp7Ny52FmjB3AwN6dc4nIvpkpXHJyDjPHDkh0KMaYTsr6HIwxpouzPgdjjDFxYcnBGGNMBEsOxhhjIlhyMMYYE8GSgzHGmAiWHOKosi7AL/61gXnv70p0KMYY0yaxrASXJiIrRWSNiGwQkXud7SNEZIWI5IvIsyISdd5pEbnL2WeLiFzYYPssZ1u+iPw4fpeUOPPe38lTH+zm169u6tCFeXYUV3Ll4+9xxv1v8eDrmwnFsD6Eab1qX4Cr/vgBJ/336zz21tZEh2NMu4il5lAHzFDVScBkYJaITAMeAB5W1VHAYeDGxgVFZBxwNTAemAU8LiJuZ1W5PwAXAeOAa5x9k9qEwVm4RMhIcTOwV8c8ABcMKVfPWc7qPaXsK63lL+/t4qkPdnXIuburV9ftZ93eUmr8QR57K5/yWn+iQzIm7ppNDhpW6bz0Ol8KzACed7bPA66IUvxyYL6q1qnqTiCf8OpxU4F8Vd2hqj5gvrNvUjt3TDbLfzKT9388k6yM5hfwiYdDVT7KavzUP8tY4w/ywY6DHXLu7mpwVvqR79NT3Em/Qp8x0cQ0fYZzp78KGEX4jn87UKqqAWeXAmBwlKKDgeUNXjfcb0+j7VEXdBWRm4CbAIYNGxZLuAnVN7NjkkK9PhleMlLc1AXCk+ile11JvyhQZzd9ZD8e/tIprNp9iC9NGUqKp/277spq/KR6XEetaGZMe4opOahqEJgsIlnAS8DY9gyq0bnnAHMgPH1GR523M1iyuYglW4oYmZ3JV04/Pur6Bx63i398cxq3PbOa4oo6LpmYw41nnZCAaLuXWRMGMmvCwA45132vbOQv7+3C7RKeuPZUZow9rkPOa7q3Fk28p6qlIrIEmA5kiYjHqT0MAfZGKbIXGNrgdcP9mtpugBdXFXD3P9dT4w+S5nWxfMchnrj2tKj7npTTi8V3nNvBEZqOUFJZx7z3dxMIKYGQcs/LGyw5mA4Ry2ilbKfGgIikA+cDm4AlwGxnt+uAl6MUXwhcLSKpIjICGA2sBD4ERjsjnlIId1ovbOO1dClz39t5ZMnQWn+If2/YT13A1lzublI8LmiwOmxmqk2kbDpGLI2lOcASEVlL+J/6YlVdBPwIuENE8oF+wFwAEblMRH4BoKobgAXARuB14BZVDTq1jVuBfxNONAucfY2jV7q34f8EPC4XHpc9ltLd9Erz8r+zJ5LdI5WR2Zn87ppTEh2S6SZsyu5OauuBCmY/8T6q4A+G+MUVE7gqd2jzBY0xppHWTNltddROasxxPXn3hzPYcqCCQVlpDOmTkeiQjDHdiCWHTqx3hpepI/omOgxjTDdkjdjGGGMiWHIwANT4gnEbDaWqHKry2RxPxiQxa1bq5jbsK+OuF9exYW85IuEpQH7zhZMZ0LN1c0OpKl+du5IPdhxkwuBevPDtM6I+vGeM6dzsr7YT2Hqggufy9rBo7T4qOnASt32lNVz1fx+wtqCMoIYfsvrP1mK+8MT7BIKhVh2zqKKO5TsOEgwpW/ZXsKOkKs5RG2M6gtUcEqikso6bnspjY2E5LhFEIBBUbvnMKP5rxihEpPmDtMG893fhCxydBAKhcJPQW5uLuHB8y6eH6N8jlZEDevDJwWr690xhWF8bZWVMMrLk0AGWbSvhf/69mb2lNUw7oR8/ufgkjuuVxlX/9wGfHKwm0Kht/ol3tpOR4uYbZ7fvHEnr95bhD0b2C9T4guQXVXLh+OaPUVhWw7f/toodJVXMHDuAB2dPYuGtZ7LtQCWjBvSwieKMSVJdvlkpEAxxyWPv8pn/fYequkDzBeJs5c5DfOOpD1lTUEZJpY9X1xVy+e/f4/X1+zlQVhuRGABnnYBt+FvZtBOrsTm98LojayfpKW5G9M+M6Ri3PbOadXvLqKgN8Pr6/Ty5bAepHjcTBve2xGBMEuvyycEXDLH1QAV7DlVTUdvxyeH3b2+j1v/pP/mQhlcS++v7O6nyNT06KKiwqbB9V5O7/ozheBt1FrtF6Jnm5bMnxTa5266D1dTnt9pAiG0HKo9dwBiTFLp8cshI8fDa7eew8NazGNi7Y1Zna2h/eW3ENl8wRM0xEgOE51pr75GgQ/tm8PdvnM6oAT1IcbvwuoXc4X144TtnxLxGwecm5pDu1BDSvC4umzyoPUM2xnSQbtHnMGpAj4Sd+8LxA9l9cMeRxXgA3C7hkomD2FlS1WTtQVUZO7Bnu8d36rA+vHnHuZRU1uF1u+id7m1R+Z9eMo7Rx/Vk8/5yzj9pIGeN7t9OkRpjOlK3SA6J9O1zR/LOliJ2FFehhNd8vvHMEdx41gie+mAXNf5gRA0h3RvujO7INvv+PVJbVc7lEq6Z2vlX6OuOqn0BXlq9l3Svm8smDbLnTUyLWHJoZ5mpHhbeehYrdh6i4HANU4b34fh+4c7e5749na/OXcmB8lp8gRAelxACZp82hNtnjk5s4CbpXTNnOVsOVCAI7+WX8NurJic6pKRT4wvyr7X72LivnKx0L5+bPIiR2YlriehIlhw6gIgw7YR+EduH9Mng7TvPZcXOQ6zZU0p6ipsLxg1MSN+ISV5l1X4efWsrr6wrxC3CF04bwjfPPoG1BWXUV0qXbClOaIzJaOnWYr7z91UoUO0L4nEJT/xnOxeOH8hDV03q8jUxSw4JVp84oiWPrkZVWb2nlPwDlQzrl8HpI/q2+4N+XV2tP8gVj7/H3sPV+JxnVuYs3cGy/BLGDerJ9qIqRIRzx2QnONLkkl9Uwbf+turIaozAkaVa39i4n18tSuHnl8fwIFASazY5iMhQ4CngOECBOar6qIhMAv4I9AB2AV9R1YixlyJyO/BNwgNw/qSqjzjbJzvl04AAcLOqrmz7JZnOqKzGz7VPrmB78adDXXN6p/HMN6cxoJfVlFrrlbWF4WbJBg8z1gVCbNlfwWNXn0JhWQ3pKR6usFFkR1FVDpTX0TczJerIvD/+Zwe+Jp4zqvWHeObDT7jjwjH0SmvZAI5kEku9KADcqarjgGnALSIyDngS+LGqngy8BPygcUERmUA4MUwFJgGXisgo5+0HgXtVdTJwj/PadFF3vbCWzfvLqfYFj3ztOljFrf/4KNGhJbXlOw9SHWXEW60/yOb95Xx1+nBmnzakyzeBtNR/v7yesx98mxm/fYdqX+TzT29tOkDwGGPJvW4XebsOtWeICdfsb4yqFqrqR873FYTXfB4MjAGWOrstBr4QpfhJwApVrXbWjf4PcGX9oYFezve9gX2tvQjTuVX7Ary5qShiqo5gCNYUlLG/LPJZEBOboX0yot75pnndHGc1sia9srYQf1A5WOkjvyjywc1YVk9u5wkMEq5FtxMiMhw4BVgBbAAud976IhBtgeP1wNki0k9EMoCLG+z3XeB/RGQP8L/AXU2c8yYRyRORvOJi61RLRlV1wXCjYhQet1BW03Ez0XY1V+UOxeOK/OF6XMLFJ+ckIKLk8J3zRuJxCScP7s1JOb0i3j/9hL4cqzvMHwxxyrCs9guwE4g5OYhID+AF4LtO38LXgZtFZBXQE/A1LqOqm4AHgDeA14GPgfo68HeA76nqUOB7wNxo51XVOaqaq6q52dnWqZaM+vdIoU9G9LZZQWKex8lEGtg7jbnXTSG7ZyoZKW7SvW6G9k3nmZumkZlq402actM5I8n/9cUs+Pb0iClkAG4+bxSpTcwSkOpxccH441r9bFCyEI2h/iQiXmAR8G9VfSjK+2OAv6vq1GaO82ugQFUfF5EyIEtVVcJDVspUNTKFN5Cbm6t5eXnNxms6n9fWFXLHgo+paTDPVLrXzd2XnMS1045PYGRdQyikbC2qwONyMTI700aBxcFzeXv46T/Xo3BkavuMFDcn5fTibzdOJSMleZKviKxS1dyWlIlltJIQvqvf1DAxiMgAVS0SERfwU8Ijj6KVr99vGOH+hmnOW/uAc4F3gBnAtpYEbpLLRSfn0Cvdy2/f2ML24iqG9Enn9pmjuaAVa0aYSC6XMHbgMe+tjgiGlPkffsL2okrOHNWfmTFOstjdfDF3KGePzubpFbtZW1BGVoaXq3KHcsbIft0i+TZbcxCRs4B3gXVA/W3fT4DRwC3O6xeBu5xawCDgSVW92Cn/LtAP8AN3qOpbDY77KOEEVUt4KOuqY8ViNQdj2u7mp1exZHMRNf4Q6V43d100lq+dMTzRYZl21JqaQ0zNSp2FJQdj2qaqLsCke984ah2RQVlpvP/jmQmMyrS31iQHG/xsTDfijjKyKdVjizKZSJYcjOlG0rxubj5vJOletzO6ycU9l45LdFimE0qe7nZjTFzcccGJnD0mm90Hq5k8NCuh652YzsuSgzHd0JThfZkyvG+iwzCdmDUrGWOMiWDJwRhjTARLDsYYYyJYcjDGGBPBkoMxxpgIlhyMMcZESKrpM0SkGNid6DiA/kBJooNoA4s/8ZL9Giz+xGvJNRyvqi1a8yCpkkNnISJ5LZ2npDOx+BMv2a/B4k+89r4Ga1YyxhgTwZKDMcaYCJYcWmdOogNoI4s/8ZL9Giz+xGvXa7A+B2OMMRGs5mCMMSaCJQdjjDERul1yEJE/i0iRiKxvsO1/RGSziKwVkZdEJKvBe3eJSL6IbBGRC5s59mMiUtngdaqIPOuUXyEiw5PwGq4XkWIR+dj5+kZnjF9E/ioiOxvEOdnZLs415TvHPjXJ4j9PRMoabL+nrfG34zWIiNwnIltFZJOI3NZgezJ8Bk3FH/fPoJ3if7dBjPtE5J8NrqvlP39V7VZfwDnAqcD6BtsuADzO9w8ADzjfjwPWAKnACGA74G7iuLnA34DKBttuBv7ofH818GwSXsP1wO87+2cA/BWYHWX7xcBrgADTgBVJFv95wKJk+DsAbgCeAlzO6wFJ9hk0FX/cP4P2iL/R8V8AvtaWn3+3qzmo6lLgUKNtb6hqwHm5HBjifH85MF9V61R1J5APTG18TBFxA/8D/LDRW5cD85zvnwdmikjkIr6d+xrirj3iP4bLgac0bDmQJSI5SRR/u2ina/gO8AtVDTnHK2pQPhk+g6bij7v2/B0SkV7ADOCfDcq3+Off7ZJDDL5OOMsCDAb2NHivwNnW2K3AQlUtbLT9SHnnQy8D+sU12ujieQ0AX3Cqo8+LyND4hhpVa+IHuM+J82ERSW1F+XiJZ/wA00VkjYi8JiLj2yHeaFpzDSOBL4lInhPr6BaWj6d4xg8d/xm09ncI4ArgLVUtb2V5wJLDUUTkbiAAPN2CMoOALwK/a6+4WqIdruFfwHBVnQgs5tOaULtoTfyOu4CxwBSgL/CjOIcWk3aI/yPC8+JMIvz5/DM+kTatDdeQCtRqeEqHPwF/jndssWiH+Dv0M2hD/PWuAZ5paxyWHBwicj1wKfAVdRrqgL1AwzvlIc62hk4BRgH5IrILyBCR/MblRcQD9AYOtkf8zjmuJ87XoKoHVbXO2e9J4LT2ib5N8aOqhU61uQ74C59Wu2MqHw/tEb+qlqtqpfP9q4BXRPq3R/xtvQbCd6QvOt+/BExsYfk2a4/4O/IzaGP8OHFNBV5psLl1P/9YOia62hcwnKM7gmYBG4HsRvuN5+iOoB003xHUsDP3Fo7ukF6QhNeQ0+D7zwPLO2P89XES7nR7BLjfeX0JR3fGrUyy+Afy6cOqU4FP6l93wmu4H/i68/15wIdJ9hk0FX+7fAbxjt/Z99vAvEbbWvXzb/MHlGxfhKtbhYCf8J3CjYQ7ePYAHztff2yw/92ERwdsAS5qsP1VYFCU4zf8x5oGPOccfyVwQhJew2+ADc4v5xJgbGeMH3gbWAesB/4O9HC2C/AHp/w6IDfJ4r+1wc9/OXBGZ/0dArII37GuAz4AJiXZZ9BU/HH/DNojfuf1O8CsRudq1c/fps8wxhgTwfocjDHGRLDkYIwxJoIlB2OMMREsORhjjIlgycEYY0wESw7GGGMiWHIwxhgT4f8BfTLP2pzNUqcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(data[:,1], data[:,0], s=erra**2)\n",
    "#plt.scatter(data[:,1], data[:,0], s=data[:,2])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 253,
   "id": "8e90c12c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 21844.45218279, -11932.64211481])"
      ]
     },
     "execution_count": 253,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.matmul(model.coef_, mytransf(np.array([(0.0,0.01,0)]))[0]) + model.intercept_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 257,
   "id": "585f4c9a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW8AAADnCAYAAADRqNcVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/SrBM8AAAACXBIWXMAAAsTAAALEwEAmpwYAAAVV0lEQVR4nO3de6xlZX3G8e/DCNKgFXAQEEalijZoy2jH8doKgjgSK9qggViLlxa12mhjY0USb20TL1Vig5WOOhEbr62iRBEYqK2aKDhQriIyIpYZkWFAxRswc/bTP/Y6uj3sc/Y6e62z93rZzydZOXtdzlq/WYHffs+73vX+ZJuIiCjLHtMOICIili/JOyKiQEneEREFSvKOiChQkndERIHuN+0AIiKm6dlH7+Pb75irdexlV919ge0NKxxSLUneETHTdt4xxyUXHFrr2D0P/t7qFQ6ntiTviJhxZs69aQexbEneETHTDPQo72XFJO+ImHk90vKOiCiKMbvSbRIRURYDc+k2iYgoT/q8IyIKY2CuwNlVk7wjYuaV1+Od5B0RM844fd4REaWxYVd5uTvJOyJmnZhD0w5i2ZK8I2KmGeil5R0RUZ60vCMiCtN/SSfJOyKiKAZ2uby6NEneETHTjJgrsKhYkndEzLye020SEVGU9HlHRBRJzKXPOyKiLP1KOkneERFFscU9XjXtMJYtyTsiZl4vfd4REWXpP7BMt0lERGHywDIiojh5YBkRUai5vKQTEVEWI3a5vFRYXsQRES3KA8uIiAIZpdskIqJEeWAZEVEYmwwVjIgoTf+BZV6Pj4goTokPLMuLOCKiRUb0XG9pQtJ7JH1H0lWSzpG0b5PzdbLlvXr/VX7Emj2nHUZEbcbTDmGi2vrXjrpvo65z8827uf2OXuOhIhNqeW8GTrO9W9K7gNOAvx/3ZJ1M3o9YsyeXXrBm2mFE1Dbn3rRDqK3XQurtMfrfO+fR1xl1nl0j7usxx9828hqjGOhN4IGl7QsHVr8JnNjkfJ1M3hERk6PllEFbLWnLwPpG2xvHuOjLgU+P8Xu/luQdETPNsJzRJjttr1tsp6SLgIOG7Drd9heqY04HdgMfX2aovyXJOyJmmq3Wuk1sH7vUfkkvBZ4LHGPX6FdaQpJ3RMy8SbykI2kD8EbgGbZ/2fR8Sd4RMdP683lPZG6TM4H7A5slAXzT9qvGPVmSd0TMuMlU0rH9qDbPl+QdETOtP1SwvFkFR37dSNokaYekawa2vU3SdklXVMvxi/zuBknXS9oq6U1tBh4R931zeMmlDfNzm9RZuqTO3wofBTYM2X6G7bXVct7CnZJWAR8AngMcAZws6YgmwUZErIQee9RaumRkNLa/CtwxxrnXA1tt32j7HuBTwAljnCciYsX0p4RVraVLmnyVvLaaYGWTpP2G7D8EuHlgfVu1LSKiUyYxMVXbxk3eHwQeCawFbgHe2zQQSadK2iJpy223zzU9XURELf1ZBfeotXTJWKNNbN86/1nSh4AvDjlsOzA4u9Sh1bbFzrkR2Aiw7si9Z2uKtoiYmv7r8d1KzHWMFbGkgwdWXwBcM+SwbwGHSzpM0l7AScC541wvImLl3Edb3pI+CRxFfzatbcBbgaMkraX/pXUT8Mrq2IcCH7Z9fDVn7WuBC4BVwCbb167EPyIiookJvWHZqpHJ2/bJQzZ/ZJFjfwgcP7B+HnCvYYQREV0xP9qkNHnDMiJmXte6ROpI8o6ImTZfw7I0Sd4RMdMM7E7LOyKiPOk2iYgoTQffnqwjyTsiZtoEizG0Ksk7ImZeWt4REYUptRhDknfECHPuTTuE2notFSiYlFEFFUbd+Tb+tUbs7uWBZUREcdLnHRFRGqfbJCKiOOnzjogo1H0yeUvaBDwX2GH7cdW29wB/CtwDfA94me2fDPndm4CfAXPAbtvrWos8IqIFRswV+MBy3Orxm4HH2f5D4LvAaUv8/tFVhfkk7ojopB6qtXTJWNXjbV9oe3e1+k36Jc4iIopjz1YB4kEvB768yD4DF0q6TNKpS50kBYgjYlps1Vq6pNEDS0mnA7uBjy9yyNNtb5f0EGCzpO9ULfl7SQHiiO7ojXw9ZjLmPCIVjNpfS/da1XWM3fKW9FL6DzJfbA+/g7a3Vz93AOcA68e9XkTESplky1vSGyRZ0uom5xm3evwG4I3A82z/cpFj9pH0wPnPwHEMrzIfETE1Nsz1VGtpStIa+rnw/5qea2TyrqrHfwN4jKRtkl4BnAk8kH5XyBWSzqqOfaik+YLDBwJfl3QlcCnwJdvnNw04IqJtExxtcgb9hm/j/p4Vqx5v+0bgyEbRRUSsMMNEHkZKOgHYbvtKqfn18oZlRMy4ZT2wXC1py8D6xmqwRf9M0kXAQUN+73TgzfS7TFqR5B0RM28Zg1Z2LvXCoe1jh22X9AfAYcB8q/tQ4HJJ623/aHnR9iV5R8TMW+luE9tXAw+ZX6+mDllne+e450zyjogVUWeseG9Ek3fU63qtFGMwRc5tkuQdETOvlXd9lnU9P6LpOZK8I2Lmde3V9zqSvCNippnuzVtSR5J3RMy8EidTSvKOiNlmcAuvvk9akndEzLx0m0REFGjSo03akOQdETNtUnObtC3JOyKWbWSRhJZMpCSEgQKTd63XiiRtkrRD0jUD2/aXtFnSDdXP/Rb53VOqY26QdEpbgUdEtMWut3RJ3XdCP8q9K8i/CbjY9uHAxdX6b5G0P/BW4En0q+i8dbEkHxExHcK9ekuX1ErewyrIAycAZ1efzwaeP+RXnw1stn2H7R8Dm7n3l0BExHS55tIhTfq8D7R9S/X5R/Qr5yx0CHDzwPq2atu9VNXlTwV42CHpio+ICXGZDyxbmUqrKkDc6HvJ9kbb62yvO+DBq9oIKyKingJb3k2S962SDgaofu4Ycsx2YM3A+qHVtoiIDlHNpTuaJO9zgfnRI6cAXxhyzAXAcZL2qx5UHldti4jojl7NpUPqDhUcVkH+ncCzJN0AHFutI2mdpA8D2L4D+AfgW9XyjmpbREQ3zI/zrrN0SK0ng4tUkAc4ZsixW4C/HFjfBGwaK7qIuE+bG9GRPDein7mtbuiujeGuI8M6IiKSvCMiCtSxLpE6krwjYuYpLe+IiMJY0LFX3+tI8o6ISMs7IqJASd4REQVK8o6IaM+uEa+ku41X1gstxpDkHREzL6NNIiJKlOQdEVGetLwjIkpUYJ/32FPCSnqMpCsGljslvX7BMUdJ+unAMW9pHHFERJvqFmJooXUu6W8kfUfStZLe3eRcY7e8bV8PrK0CWkW/yMI5Qw79mu3njnudiIgVN4FuE0lH06/9e6TtuyU9pMn52uo2OQb4nu0ftHS+iIiJ0WQKLbwaeKftuwFsD6s+VlsrNSyBk4BPLrLvKZKulPRlSY9d7ASSTpW0RdKW226faymsiIga6nebrJ7PU9Vy6jKu8mjgjyVdIul/JD2xSciNW96S9gKeB5w2ZPflwMNt/1zS8cDngcOHncf2RmAjwLoj9y7w2W+Uas4dq2+1hF5BY9pGFVqA0ZXFdnnp9mUbd0Ne1miTnbbXLXou6SLgoCG7Tqefb/cHngw8EfiMpN+rCrgvWxvdJs8BLrd968Idtu8c+HyepH+VtNr2zhauGxHRjpZGm9g+drF9kl4NfK5K1pdK6gGrgdvGuVYb3SYns0iXiaSDJKn6vL663u0tXDMioj2TGW3yeeBoAEmPBvYCxm7INmp5S9oHeBbwyoFtrwKwfRZwIvBqSbuBXwEnjfsnQkTESpnQSzqbgE2SrgHuAU5pkg8bJW/bvwAevGDbWQOfzwTObHKNiIgV5cmMNrF9D/DnbZ0vb1hGRBTYH5DkHRGR5B0RUZ5MTBURndcbObp6cuZGPK+7y6uW3N9roxhDoZK8IyLS8o6IKMyERpu0Lck7IiIt74iIsog8sIyIKFOSd0REYZY3q2BnJHlHROSBZUREedLyjoiZUOdFn16NCfNG1cz6hfcaEUdLL+kUmLwbz+ct6SZJV1fV4bcM2S9J/yJpq6SrJD2h6TUjIlozwerxbWqr5X30EtVxnkO/9NnhwJOAD1Y/IyI6ocRuk7YKEC/lBOBj7vsmsK+kgydw3YiIegpsebeRvA1cKOmyRSopHwLcPLC+rdr2W1I9PiKmRb16S5e00W3ydNvbJT0E2CzpO7a/utyTpHp8RExFB1vVdTRuedveXv3cAZwDrF9wyHZgzcD6odW2iIip0zKWLmmUvCXtI+mB85+B44BrFhx2LvAX1aiTJwM/tX1Lk+tGRLSqwD7vpt0mBwLnSJo/1ydsn7+ggvx5wPHAVuCXwMsaXjMiolUljjZpWj3+RuDIIdsHK8gbeE2T60TEZI2qcNOWUc8Af9bbe+nf9+y+pJM3LCNitqUYQ0REodLyjogoz8z1eUdE3CckeUdElCct74iI0pgUY4iIKE0KEEdELNPciKT5o137Lrl/l1tKYRNI3pLWAmcBewO7gb+2fem455vElLAREZ0mu9bS0LuBt9teC7ylWh9bWt4RMdsmN2+Jgd+tPj8I+GGTkyV5R8TMW0af9+oF5R43VtNZ1/F64AJJ/0y/1+Opta86RJJ3RMy8Zbwev9P2ukXPI10EHDRk1+nAMcDf2v6spBcBHwGOXWaovzZ28pa0BvgY/ZkFTf8b6P0LjjkK+ALw/WrT52y/Y9xrRkSsiJa6TWwvmowlfQx4XbX6H8CHm1yrSct7N/AG25dXc3pfJmmz7W8vOO5rtp/b4DoRESvHExsq+EPgGcB/A88EbmhysrGTd1VQ4Zbq888kXUe/NuXC5B0R0W2TSd5/Bbxf0v2Au4BhNX9ra6XPW9IjgMcDlwzZ/RRJV9L/1vk729e2cc2IiDZM6iUd218H/qit8zVO3pIeAHwWeL3tOxfsvhx4uO2fSzoe+Dxw+CLnOZXqm+hhh+Q5akTp5mo0Z3eNqAz53buGPfv7jbu857JiWox65b1i2bSG5Z70E/fHbX9u4X7bd9r+efX5PGBPSauHncv2RtvrbK874MGrmoQVEVFf3fqVHcvvTUabiP5Ql+tsv2+RYw4CbrVtSevpf1ncPu41IyJWwqxV0nka8BLgaklXVNveDDwMfl3H8kTg1ZJ2A78CTqpqWkZEdEeBWanJaJOvw9IdVrbPBM4c9xoREZOQWQUjIkpjoMAOgSTviJh5s9bnHRFRvBRjiIgokZ1uk4iumXM5fw/3ChryUOcFnDp3fpeXftVk882PWXL/nfd8vcZVRkvLOyKiREneERHlScs7IqI0ZnQl5A5K8o6ImZeWd0REiTLaJCKiPGl5R0SUpoPTvdaR5B0RK2KuRlfEqHHeBzzv+iX33+i7lxXTMAJU4APLpsUYNki6XtJWSW8asv/+kj5d7b+kKpcWEdEpsmstXTJ28pa0CvgA8BzgCOBkSUcsOOwVwI9tPwo4A3jXuNeLiFgRhVbSadLyXg9stX2j7XuATwEnLDjmBODs6vN/AsdUFXgiIjrCv5nfZNTSIU2S9yHAzQPr26ptQ4+xvRv4KfDgYSeTdKqkLZK23Hb7XIOwIiKWR663dEmjPu82pQBxRExNgS3vJqNNtgNrBtYPrbYNO2abpPsBDyIFiCOiSzx7o02+BRwu6TBJewEnAecuOOZc4JTq84nAf6UAcUR0ToEPLJsUIN4t6bXABcAqYJPtayW9A9hi+1zgI8C/S9oK3EE/wUdEdErXhgHW0eglHdvnAect2PaWgc93AS9c7nkvu+runasO3vqDgU2rgZ3jxjlhiXVllBQrlBVvh2O9eeGGhbE+vJXLzFryXim2Dxhcl7TF9rppxbMciXVllBQrlBXvzMdq6pX96ZhOJu+IiEkR3Xt7so7ODBWMiJiaXq/e0oCkF0q6VlJP0roF+06rphG5XtKz65yvlJb3xmkHsAyJdWWUFCuUFe9sxzq5bpNrgD8D/m1wYzWtyEnAY4GHAhdJerTtJd9WLCJ52y7mP67EujJKihXKijexTma0ie3rAIbMEHIC8CnbdwPfr0bnrQe+sdT50m0SEVH/DcvV89N4VMupLVy9zlQj99Lp5D1qytmukXSTpKslXSFpy7TjGSRpk6Qdkq4Z2La/pM2Sbqh+7jfNGOctEuvbJG2v7u0Vko6fZozzJK2R9BVJ3676M19Xbe/cvV0i1q7e270lXSrpyiret1fbD6ummN5aTTm9V7MrLWtiqp3z03hUy2/9JSDpIknXDFkWTtrXWGeTd80pZ7voaNtrOzj06qPAhgXb3gRcbPtw4OJqvQs+yr1jBTijurdrq3cMumA38AbbRwBPBl5T/XfaxXu7WKzQzXt7N/BM20cCa4ENkp5Mf2rpM6qppn9Mf+rp8c1Xj6+zjDqVfaztxw1ZvrDEr9WZauReOpu8qTflbNRk+6v033IdNDhl79nA8ycZ02IWibWTbN9i+/Lq88+A6+j/ydu5e7tErJ3kvp9Xq3tWi4Fn0p9iGlq6t1MuxnAucFJVvOYw4HDg0lG/1OXkPVY/0JQZuFDSZS31ha20A23fUn3+EXDgNIOp4bWSrqq6VabeDbFQVSnq8cAldPzeLogVOnpvJa2SdAWwA9gMfA/4STXFNLSVFyYwq6CkF0jaBjwF+JKkC/qX9rXAZ4BvA+cDrxk10gS6nbxL9HTbT6Df1fMaSX8y7YDqqiYM6/KbCh8EHkn/z+dbgPdONZoFJD0A+Czwett3Du7r2r0dEmtn763tOdtr6XclrAd+v/2LAD3XW5pcxj7H9qG272/7QNvPHtj3T7Yfafsxtr9c53xdTt5j9QNNk+3t1c8dwDn0/2PrslslHQxQ/dwx5XgWZfvW6n/kHvAhOnRvJe1JPxl+3Pbnqs2dvLfDYu3yvZ1n+yfAV+i3WvetppiGVvLC7FXSWWl1ppztDEn7SHrg/GfgOPqD8rtscMreU4ClHqpM1XwirLyAjtxb9QftfgS4zvb7BnZ17t4uFmuH7+0BkvatPv8O8Cz6/fRfoT/FNLR1bwtM3p19SWexKWenHNZSDgTOqQbg3w/4hO3zpxvSb0j6JHAU/XGq24C3Au8EPiPpFcAPgBdNL8LfWCTWoyStpf9H7k3AK6cV3wJPA14CXF31zQK8mW7e28ViPbmj9/Zg4Oxq5NkewGdsf1HSt4FPSfpH4H/pfyGNz8BceTNTKbURImKWPej+B/qpD31xrWPPv+mMy7oyDLizLe+IiIkpsBGb5B0Rs21+tElhkrwjItLyjogoUJJ3RERhbJgb+UJj5yR5R0Sk5R0RUaAk74iI0jSft2QakrwjYrYZ+tO6lCXJOyKiwNfjk7wjYrbZ0EvyjogoTx5YRkSUx2l5R0SUpntzddeR5B0Rsy0TU0VElMeA83p8RERhbMg474iI8jjdJhERBSqw5Z0alhEx0ySdD6yuefhO2xtWMp66krwjIgq0x7QDiIiI5UvyjogoUJJ3RESBkrwjIgqU5B0RUaD/B82rZnaquDFMAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "lat=(np.arange(29.96, 30.13, 0.01)-30.05)\n",
    "lon=(np.arange(120.38, 120.69, 0.01)-120.5)\n",
    "L, B = np.meshgrid(lon, lat)\n",
    "BL0 = np.array([B, L, np.zeros_like(B)]).transpose((1,2,0)).reshape((-1,3))\n",
    "BL1 = np.array([B, L, np.ones_like(B)]).transpose((1,2,0)).reshape((-1,3))\n",
    "a = model.predict(mytransf(BL0)).reshape((len(lat), len(lon), 2))\n",
    "b = model.predict(mytransf(BL1)).reshape((len(lat), len(lon), 2))\n",
    "plt.imshow(np.diff(a[:,:,0], axis=1)/a[:,:-1,0], origin='lower')\n",
    "plt.colorbar()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 296,
   "id": "cae65eef",
   "metadata": {},
   "outputs": [],
   "source": [
    "def parttofull(x):\n",
    "    i, p, l, p2, lp, l2, p3, p2l, l2p, l3, h = x.transpose()\n",
    "    lh = ph = h2 = plh = h2l = h2p = l2h = p2h = h3 = np.zeros_like(i)\n",
    "    return np.array([i, l, p, h, lp, lh, ph, l2, p2, h2, plh, l3, p2l, h2l, l2p, p3, h2p, l2h, p2h, h3])\n",
    "\n",
    "def write_rpc(fn, model):\n",
    "    with open(fn, 'w') as f:\n",
    "        f.write('ERR_BIAS: 0.0\\n'\n",
    "                'ERR_RAND: 0.0\\n'\n",
    "                f'LINE_OFF: {-model.intercept_[1]}\\n'\n",
    "                f'SAMP_OFF: {model.intercept_[0]}\\n'\n",
    "                'LAT_OFF: 30.05\\n'\n",
    "                'LONG_OFF: 120.5\\n'\n",
    "                'HEIGHT_OFF: 1\\n'\n",
    "                'LINE_SCALE: 1\\n'\n",
    "                'SAMP_SCALE: 1\\n'\n",
    "                'LAT_SCALE: 1\\n'\n",
    "                'LONG_SCALE: 1\\n'\n",
    "                'HEIGHT_SCALE: 1\\n')\n",
    "        line_num_coeff = -parttofull(model.coef_[1])\n",
    "        line_den_coeff = np.zeros(20)\n",
    "        line_den_coeff[0] = 1\n",
    "        samp_num_coeff = parttofull(model.coef_[0])\n",
    "        samp_den_coeff = np.zeros(20)\n",
    "        samp_den_coeff[0] = 1\n",
    "        for i in range(20):\n",
    "            f.write(f'LINE_NUM_COEFF_{i+1}: {line_num_coeff[i]}\\n')\n",
    "        for i in range(20):\n",
    "            f.write(f'LINE_DEN_COEFF_{i+1}: {line_den_coeff[i]}\\n')\n",
    "        for i in range(20):\n",
    "            f.write(f'SAMP_NUM_COEFF_{i+1}: {samp_num_coeff[i]}\\n')\n",
    "        for i in range(20):\n",
    "            f.write(f'SAMP_DEN_COEFF_{i+1}: {samp_den_coeff[i]}\\n')\n",
    "write_rpc('/mnt/2thdd-btrfs/rs/USGS_Declass/D3C1215-401419A011/D3C1215-401419A011_e_RPC.TXT', model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a8f94fa6",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
